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Abstract. We study statistical properties of a zero-range process (ZRP) on random networks. We derive 
an analytic expression for the distribution of particles (also called node occupation distribution) in the 
steady state of the ZRP in the ensemble of uncorrelated random graphs. We analyze the dependence of this 
distribution on the node-degree distribution. In particular, we show that when the degree distribution is 
tuned properly, one can obtain scale-free fluctuations in the distribution of particles. Such fluctuations lead 
to a power law in the distribution of particles, just like in the ZRP with the hopping rate u(m) = 1 + b/m 
on homogeneous graphs. 
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1 Introduction 

Many statistical systems are defined on random networks 
or random lattices. Usually, this means that one considers 
a system on a lattice with some irregularities, or even on 
a purely random network, but instead of looking at what 
happens on a single network, one performs an annealed av- 
erage over a statistical ensemble of networks that is over a 
collection of random graphs with some statistical weights. 

The question of how averaging over the disorder influ- 
ences the statistical properties of the system has been pre- 
viously addressed in the context of two-dimensional sta- 
tistical models. It was shown that additional degrees of 
freedom related to fluctuations of the geometry can lead 
to quite distinct behavior, in comparison to analogous sys- 
tems defined on a fixed lattice. For example, critical prop- 
erties of the two-dimensional Ising model on a fixed lattice 
are described by Onsager exponents pQ while on a ran- 
dom lattice, represented as a sum over planar networks, 
by "dressed" KPZ-DDK ones |2p] . A similar change of 
critical exponents is observed for other models. 

In this paper we address the same problem but for 
a statistical system on complex random networks. As a 
particularly simple but interesting example we consider 
a zero-range process (ZRP), which has been thoroughly 
studied on fixed networks [4,5,6,7,8,9,10]. The particu- 
lar feature of this out-of-equilibrium model is that under 
certain conditions particles tend to condense on a single 
node. Here we analyze the influence of the averaging over 
random networks on the distribution of particles in the 
steady state. We shall discuss an ensemble of networks 



with given probability distribution of degrees. To make 
things as simple as possible we shall restrict the discussion 
to uncorrelated networks with independent node degrees. 

The paper is organized as follows. In the next section 
we shall recall the definition of the ZRP on a network and 
introduce quantities describing the steady state of this 
model. Then we shall consider the free ZRP process being 
a particular class of ZRP, in which the hopping rates do 
not depend on the distribution of particles but only on the 
connectivity of the network. For this case we will present in 
Sees. 3 and 4 an exact solution to the problem of how the 
distribution of particles averaged over nodes depends on 
the node-degree distribution and how to choose the latter 
one to obtain a scale-free distribution of particles. In the 
following Sec. 5 we will discuss finite-size effects observed 
in the distribution of particles for some models of random 
networks like Erdos-Renyi graphs. We will finish in Sec. 6 
with some concluding remarks. 



2 Model description 

Consider the ZRP on a connected simple graph with TV 
nodes and a sequence k — {fci, . . . , k^} of node degrees. 
The state of the ZRP is given by the distribution of par- 
ticles m = {mi, . . . , mpf} on nodes of the network, where 
rrii > is the number of particles at node i. The total 
number of particles M = m\ + . . . + is conserved dur- 
ing the process. A particle can hop from node i to one 
of its neighbors with rate u(rrii)/ki. The function u(m), 
called hopping rate, depends only on the number of parti- 
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cles at the departure node i. The factor takes care of 
distributing the outflow of particles equally between the 
neighbors. The hopping rates u(m) are non-negative and 
identical for all nodes. 

The ZRP is known to have a unique steady state [TT] . 
Static properties of this state are described by a parti- 
tion function Z(N, M, k) depending only on the degree 
sequence k: 



M 



N 



Z(N,M,k) = 3 

J[p{rrii)kp, 



mi,.. . ,mjv — 



with statistical weights p(m) defined as 



m 

^H = n^ Ko) = i. 



(i) 



(2) 



We will refer to Z(N, M, k) as to a microcanonical parti- 
tion function. The main quantity describing the system is 
the distribution of particles 7r(m), averaged over all config- 
urations m with the weight given by the partition function 
([TJ and over all nodes: 



(3) 



where the argument fe means that it is calculated for a 
single network with a given sequence of degrees. It is also 
called node occupation distribution. It can be calculated 
[TO] as follows: 

where Zi is the partition function for a graph with N — 
1 nodes and degrees k' = {k\, . . . , . . . , fc/v}. 

Equation (U) holds for any connected graph with a given 
degree sequence. 

Suppose now that we are interested in the behavior 
of the ZRP on a random network. In this case we have 
to take the average over networks. Denote by P{k) = 
P(ki, . . . , fc/v) the probability of choosing a network with 
the degree sequence k±, . . . , kjq- We can now define a canon- 
ical partition function as the average over all degree se- 
quences: 



Z(N, M) = p (. k ) z ( N > M > k ) 



(5) 



In general, P(k) may have a complicated form. We shall 
restrict our attention to uncorrelated networks [T21IT3] for 
which P(k) is a product measure: P{k) = il(fci) • • ■ 77(fc/v). 
This means that node degrees are independent of each 
other and that the observed degree distribution is 7T(fc)13 



1 We neglect the fact that the total number of links is often 
fixed, which leads to an additional Kroenecker delta constraint 
$2L,k 1 + ...+k N - This constraint can be usually neglected in the 
thermodynamic limit. 



The canonical partition function assumes then a simple, 
symmetric form: 



N 



Z(N,M)= Y 8 mi +...+m N ,M Hftrm), (6) 



(.= 1 



where 



pirn) 



P (m)J2n(k)k r 



(7) 



fc=i 



is an effective weight for a node occupied by m particles. 
As we see, the effective weight p(m) is calculated from the 
node degree distribution II{k) and the occupation weight 
p(m). The effective partition function Z(N, M) in Eq. ([6]) 
has the form of a partition function of the balls-in-boxes 
model with identical weights, which has been thoroughly 
studied [TH[T5, 16J. Z{N,M) is invariant with respect to 
any permutation a of node occupation numbers m, — * 



3 Free ZRP 

In this section we shall consider a particular example of a 
ZRP for which the hopping rate u(m) = 1 is independent 
of m. We shall call it free ZRP (FZRP). In this case, also 
the occupation weight ([2]) is constant, p(m) = I, and the 
canonical partition function ([6]) reduces to 

TV 

i+...+mjv,M HKmi), (8) 



mi,...,mjv 



where 



(9) 



k=l 



is the m-th moment of the node-degree distribution. The 
probability that a node is occupied by m particles is now: 

. , Z(N-l,M-m) , s 

TT m = — i — ■ -u(m), (10) 

V ' Z(N,M) ^ y h y ' 

just as it was in the balls-in-boxes model [14]. 

For further convenience, let us introduce a generating 
function for the moments fj,(m): 



z 



(11) 



which encodes the same information as il(fc). Indeed, in- 
serting Eq. ([9|) into the last equation we see that it can 
be interpreted as a Fourier transform of the node-degree 
distribution, 



M(-iz) = ^2n(k) e - tkz . 



(12) 



From the generating function one can formally reconstruct 
the moments, 



MM = ir- <t> M(z)z- m ~ 1 dz, 

ZTTl ' 



(13) 
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as well as the degree distribution, 

n(k) = — [ dze izk M(-iz) 

271 J-7, 



(14) 



The partition function © is well defined if all moments of 
the distribution II(k) are finite. Usually, we are interested 
in the behavior of the system in the thermodynamic limit 
N — > oo. We can distinguish two cases: (a) the limiting dis- 
tribution II(k) for N — > oo has all moments finite, as for 
instance for Erdos-Renyi (ER) graphs, where it is Poisso- 
nian, (b) higher moments of the limiting distribution ZZ (k) 
diverge for N — > oo as it happens for scale- free graphs [17] . 
For (a), the large ./V limit presents no difficulty, while for 
(b) it has to be taken very carefully since it depends on 
the details of how the ensemble is defined. Moreover, in 
case (a) one can show that for random graphs^] the prob- 
ability of any sequence of degrees k factorizes in the limit 
N — > co [17]. This factorization often breaks down for 
(b). One observes particularly strong deviations from the 
factorization for II(k) ~ k 1 with 2 < 7 < 3 where finite- 
size effects are especially strong [T2l[T51[T^] . Below we shall 
discuss only the case (a) which is free of these problems. 



4 Power-law distribution of particles 

We have shown in the previous section that averaging over 
fluctuating geometries leads to an effective model with the 
partition function ([5]) and weights /x(m) being the mo- 
ments of the node-degree distribution. This model has an 
interesting critical behavior for weights which fall off like 



/i(m) ~ m . For example [20], for 

/i(m) oc r(m + l)/r(m + b + 1) ~ m 



-b 



(15) 



one observes a condensation of particles when the density 
of particles p = M/N is larger than a critical density p c = 
1/(6 — 2). In the thermodynamic limit, at the critical point 
p = pc fluctuations of the number of particles become 
scale-free and w(m) = /x(m) displays a power law. Below 
p c it has an exponential cut-off: 



7r(m) = p(m) exp(A — Bm), 



(16) 



where the constants A and B are chosen so that 
the normalization Y] m 7r(m) = 1 and the density of parti- 
cles ^2 m to7t(to) = p are fixed. Above p c , the distribution 
7r(m) is approximately given by p(m) but with an addi- 
tional peak centered around to* rs M — p c N. 

One now can ask whether the weights (TT5|) can be ob- 
tained in our FZRP by tuning the node-degree distribu- 
tion of the underlying network. Before we proceed, it is 
important to notice that the model given by the partition 
function ([5]) is invariant with respect to the rescaling: 



(17) 



By a random graph we understand here a graph being max- 
imally random among all graphs with a given sequence of de- 
grees. 



Indeed, the partition function ([8]) changes only by a factor: 
Z(N, M) -> Af N (/) M Z{N, M), which is constant for given 
N and M, while physical quantities stay intact because the 
normalization factor cancels out. Thus, we expect that if 
the moments are given by 



. . r(m+l) ,_ 

1 (771+1 + 0) 



(18) 



then the degree distribution of node occupation numbers 
at the critical density should be given by 



7r(m) oc 



r(m + l) 
r(m + 6 + 1) 



(19) 



The question we face now is whether there is a node- 
degree distribution n(k) which has moments given by 
Eq. (|18p. First of all, we observe that the parameter <f> in 
Eq. (fT8|) plays the role of a scale parameter of the distri- 
bution n{k) as follows from the definition of the moments 
([9]): under the rescaling k — > k/<j> the moments change as 
p{m) — > p(m)4> rn . We will use the freedom of choosing 
the parameter <fi to fix the average degree k and thus also 
the number of links L = kN/2. The parameter Af has 
to be chosen in such a way that n(k) takes the proper 
normalization of a probability. 

The integral in Eq. ([T4|) is hard to calculate and can- 
not be easily expressed in terms of elementary functions. 
However, if we assume ^ > 1, then the function M{— iz) 
goes to zero sufficiently fast when z — > ±00 and thus we 
can extend the limits of integration to ±00. In this case 
the integral can be done analytically. Equation ([14]) be- 
comes a Fourier transform of the function M(— iz) which 
is a special case of the Mittag-Leffler functions having a 
known form of an infinite series expansion (see e.g. |21|). 
Changing variables k — > x<p we obtain 



n(x<i>) = 



(-izy 



^ r{ m + 1 + b) ' 

m— v ' 



(20) 



According to Eq. (B21) of Ref. [21] , the above integral 
yields 



d> ' m\r(b — to) 

T m=0 v ; 



and hence 



Af 



m\r{b- 



TO 



(<f>-k) 



6-1 



r{b)<p b 



(21) 



(22) 



In Fig. [T] we compare n(k) computed numerically using 
the integral (fT4]) with the original limit ±7r, and calcu- 
lated by means of Eq. ([2"2"]l . Because the probability 7T(fc) 
must be non-negative, the above solution is physical only 
for k < 4> and we have to set n(k) — for k > cf>. We see 
that the integer part of cj> can be interpreted as the maxi- 
mal degree which can exist in the network. The existence 
of the upper cut-off in the node-degree distribution is not 
only a property of the approximate solution. Also when 
one uses the exact relation (fT4"]) to calculate the degree 
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k is almost linear: 



Fig. 1. n(k) calculated from the exact formula (|14[1 (points) 
and approximated one (Eq. (|22p . thick line), for TV = 2ir and 
6 = 3. Squares: = 6, circles: <fi = 30. The approximate solu- 
tion diverges for x — k/cj> > 1 and has to be cut. For < x < 1 
the approximation is the better, the larger is 0. 



a, 10 




1000 



Fig. 2. Desired (dashed line) versus real distribution of par- 
ticles n(m) for networks with finite average degree k (solid 
lines), obtained from Eqs. (|22[1 and (|9} for 6 = 3. Lines 
from left to right: = 5,10,20,40,80 which corresponds to 
k = 1.67,2.89,5.4,10.4,20.4 from Eq. JUJ. These plots ap- 
proximate also Tv(m) at the critical point. The parameter 
grows almost linearly with k. 



distribution 77 (k) for the moments of the form (fT5|) . one 
obtains negative values of II(k) for k > cj), so again one 
has to cut off the solution and set 77(fc) to zero for k > 0. 
If one now calculates moments for the distribution (1221) 
with the cut-off directly from the definition ((9]), one will 
see that they deviate slightly from those given by Eq. (jTHJ) • 
However, the deviation decreases when increases and fi- 
nally becomes negligible for sufficiently large 0. In Fig. [2] 
we plot the moments /z(m) from Eq. (Tig)) and those cal- 
culated from Eq. ([9]), for various 0. As increases, the 
curves tend asymptotically to a power law. 

The parameter is related to the average degree as 
k — Y^k=i n(k)k. For large 0, the relation between and 



k = 



\b-l 



k !}{<t>-k)k 



6-1 



dk 



l£=i(*- 



\b-i 



6+1 
.(23) 

This implies that in order to obtain the power-law distri- 
bution (fT9|) for N — > oo, the value of k should increase 
to infinity. For sparse networks k would be finite. We see 
thus that the price to pay for having a scale-free distri- 
bution of particles is to make networks denser when their 
size increases. 

The normalization factor M must be chosen 
so that the degree distribution is normalized to unity: 
X)fe=i n(k) = 1. For example, for 6 — 2,3,4 we obtain 
the following degree distributions n(k) for < k 



< 



b = 2 
b = 3 
6 = 4 



2(0 -k) 



-1)' 



■1)(20-1)' 
4(0- fc) 3 



l) 2 ' 



(24) 
(25) 
(26) 



and zero for both k = and k > <p, with given by the 
following formulas for b = 2,3: 



3k- 1, 



= ( -1 + 4fc + \J\ - 16fc+ 16fc 2 ) /2, 



and by the solution of a cubic equation for 6 = 4: 
(l + 0)(30 2 -2) 



150(0-1) 



= k. 



(27) 
(28) 



(29) 



How does it come about that the power laws are observed 
in the distribution of particles 7r(m) when one averages it 
over networks while they are not seen in 7r(m)'s for indi- 
vidual nodes, for any single network in the ensemble? The 
answer is that the effective distribution n(m) averaged 
over networks is a subtle result of a well-tuned superposi- 
tion of many exponential distributions: for a node with de- 
gree k, the distribution of particles is n(m) oc (fc//c max ) m , 
where fc max is the maximal degree in the network [9]. On 
the node with maximal degree, however, there is a conden- 
sation just as for scale-free networks [5], but it disappears 
in the thermodynamic limit. This happens because the 
critical density for the condensation becomes larger than 
p = 1/(6 — 2) which we assumed to hold in our system, 
and the system is always in the fluid phase. 



5 Other random graphs 

In the previous section we found a node degree distribu- 
tion for the ensemble of random, uncorrelated networks for 
which the corresponding FZRP has a power-law particle 
distribution n(m). What happens with FZRP on generic 
random networks? Can the particle distribution be scale 
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Fig. 3. The particle distribution 7r(m) for ER graphs, N = 
400, obtained from Eq. © by using a Poissonian degree dis- 
tribution with fc max ~ 22, calculated for n = 200 samples as 
described in the main text. The almost straight line on the 
log-log plot explains partially the quasi-power-law observed in 
numerical experiments (see the data points). 

free? Let us begin with what happens in the limit N — > oo. 
Consider some typical examples of graphs like random 
trees [22] or ER graphs for which the limiting shapes of 
the degree distribution are known: 7T(fc) = e~ 1 /(fc — 1)! 
and n(k) = e~ k k k /k\, respectively. So we can calculate 
the corresponding critical distribution 7r(m). In the first 
case, the generating function has a closed form M{z) = 
exp(z + e z ), as follows from Eq. (fTT|) . and we can de- 
duce the coefficients fi(m) from the inverse Laplace trans- 
form (fT3|) . Using the saddle-point method and integrating 
around zq log(m/ logm) we can find the asymptotic 
behavior for large m: 

log/x(m) = m(logm — log logm) + 0{m). (30) 

We see that /i(m) grows over-exponentially for large m. 
This means that in the thermodynamic limit the conden- 
sation always happens, regardless of the density of parti- 
cles. The distribution 7r(m) in the bulk falls faster than 
any power law. Similarly, one can estimate that for ran- 
dom ER graphs M(z) oc exp(fce z ) — 1 and the leading term 
in log/i(m) is also m log to. So again it is clear that in the 
limit M — » oo one cannot obtain a power-law distribution 
of particles. 

Surprisingly, in Monte Carlo (MC) simulations of ran- 
dom trees and ER graphs one observes distributions ir(m) 
which very much resemble those with power-law tails. Why 
does it happen? The answer is that this is caused by finite- 
size effects. When one repeats the calculations for finite 
systems, including sub-leading terms and a finite-size cut- 
off, one obtains a line like that in Fig. [3l which very much 
resembles a scale-free plot. The line is compared to MC 
simulations. In our analytic calculations aimed at mim- 
icking the result of the simulation, finite-size effects were 
taken into account as follows. In order to compute the mo- 
ments we introduced a finite-size cut-off to II(k). The po- 
sition of the cut-off fc max was estimated from the condition: 



JI(fc max ) w l/(nN), where n is the number of samples in 
MC simulations. Next, for the distribution with the cut-off 
fcmax we calculated the moments (|9]) and we got rid of the 
leading exponential behavior of /Lt(m), using the freedom 
(|17p . by multiplying them by a factor exp(— mB) with B 
appropriately chosen. As we can see in Fig. [3] we obtained 
almost a straight line in the log- log plot. The discussion 
from the beginning of this section tells us, however, that 
this quasi-power-law behavior is only a finite-size effect 
which will disappear in the limit N — » oo. 



6 Conclusion 

We analytically investigated the influence of annealed av- 
eraging over random networks on the statistical proper- 
ties of ZRP. In particular, we calculated how the parti- 
cle distribution depends on the node-degree distribution. 
We showed that by tuning the node-degree distribution 
we can make the corresponding zero-range process critical 
and having a power-law distribution of particles. 

We believe it is the first step towards the analysis of 
more complex systems, where the topology and the dy- 
namics of the system are coupled to each other and in- 
fluence mutually. An example of this type of interactions 
was discussed in the context of 2d statistical systems on 
random lattices (2d gravity), where a back-reaction of the 
system on the lattice was observed [2"ll3"ll2"3"] . which mani- 
fested as a change of fractal properties of the underlying 
geometry when the system became critical. It would be 
very interesting to see such an adaptation mechanism also 
for random complex networks. 
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